4 Making Maps
“Maps invest information with meaning by translating it into visual form.” — Susan Schulten
How the BBC uses R. https://medium.com/bbc-visual-and-data-journalism/how-the-bbc-visual-and-data-journalism-team-works-with-graphics-in-r-ed0b35693535
4.1 Maps with tmap
So far we’ve used the plot() method as a first-order check on our geocomputations. It is possible to make a publishable quality map using plot() but there is a lot of trial and error in getting it to the quality needed.
An efficient alternative is to use functions from the tmap package. There are other packages for making nice maps in R with some of them listed in the syllabus but I like the tmap package because it integrates sf, sp, and raster objects seamlessly.
Like ggplot2, tmap uses the ‘grammar of graphics’. The grammar separates the input data frame from the aesthetics (how data are visualised). The functions translate the data into aesthetics. The aesthetics can include the location on a geographic map (defined by the geometry), color, and other visual components.
All tmap maps start with the tm_shape() function that takes as input a spatial data frame. The function is followed by one or more layers such as tm_fill() and tm_dots() that defines how a property in the data gets translated to an aesthetic.
Consider the New Zealand simple feature data frame (nz) from the spData package. Make the data frame available, check its class, and make a plot of the geometry.
## [1] "sf" "data.frame"
The geometry column is labeled geom.
Check the native coordinate reference system (CRS) of the spatial data frame with the st_crs() function from the sf package.
## Coordinate Reference System:
## User input: EPSG:2193
## wkt:
## PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",
## GEOGCS["NZGD2000",
## DATUM["New_Zealand_Geodetic_Datum_2000",
## SPHEROID["GRS 1980",6378137,298.257222101,
## AUTHORITY["EPSG","7019"]],
## TOWGS84[0,0,0,0,0,0,0],
## AUTHORITY["EPSG","6167"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4167"]],
## PROJECTION["Transverse_Mercator"],
## PARAMETER["latitude_of_origin",0],
## PARAMETER["central_meridian",173],
## PARAMETER["scale_factor",0.9996],
## PARAMETER["false_easting",1600000],
## PARAMETER["false_northing",10000000],
## UNIT["metre",1,
## AUTHORITY["EPSG","9001"]],
## AUTHORITY["EPSG","2193"]]The CRS is called the New Zealand transverse Mercator. The distance unit is meter (https://epsg.io/2193).
To make a tmap map we first identify the spatial data frame with the tm_shape() function and then add a borders layer.
The borders separate New Zealand into 16 administrative regions.
The function tm_shape() and its subsequent drawing layers (here tm_borders()) as a ‘group’. The data in the tm_shape() function must be a spatial object of class simple feature, raster, or an S4 class spatial object.
Here we use a fill layer instead of the borders layer.
Here we layer using the fill aesthetic and then add a border aesthetic.
Layers are added with the + operator and are functionally equivalent to an overlay.
We can assign the resulting map to an object. For example here we assign the map of New Zealand to the object map_nz.
## [1] "tmap"The resulting object is of class tmap.
New spatial data are added with + tm_shape(new_object). In this case new_obj represents a new spatial data frame to be plotted over the preceding layers. When a new spatial data frame is added in this way, all subsequent aesthetic functions refer to it, until another spatial data frame is added.
For example, let’s add an elevation layer to the New Zealand map. The elevation raster (nz_elev) spatial data frame is in the spDataLarge package on GitHub.
The install_github() function from the devtools package is used to install a package on GitHub. GitHub is a company that provides hosting for software development version control using Git. Git is a distributed version-control system for tracking changes in code during software development.
Note, I’ve done this for you below (Do not run the code chunk below).
Make the data available.
Next identify the spatial data for the the new layer by adding tm_shape(map_nz). Then add the raster layer with the tm_raster() function and set the transparency level to 70% (alpha = .7).
The new map object map_nz1 builds on the existing map object map_nz by adding the raster layer nz_elev representing elevation.
We can create new layers with functions. For instance, a function like st_union() operates on the geometry column of a simple feature data frame.
As an example, here we create a line string layer as a simple feature object using three geocomputation functions. We start by creating a union over all polygons (regions) with the st_union() function applied to the nz simple feature object. The result is a multipolygon defining the coastlines.
Then we buffer this multipolgyon out to a distance of 22.2 km using the st_buffer() function. The result is a single polygon defining the coastal boundary around the entire country.
Finally we change the polygon geometry to a line string geometry with the st_cast() function. The default in st_cast() is to simplify the geometry (e.g., a LINESTRING is simpler than a POLYGON).
The operations can be linked together with the pipe (%>%) from the dplyr package.
library(dplyr)
nz_water <- st_union(nz) %>%
st_buffer(22200) %>%
st_cast(to = "LINESTRING")
class(nz_water)## [1] "sfc_LINESTRING" "sfc"Now we add the resulting spatial object as a layer to our map.
Finally, lets create a layer representing the country elevation high points (stored in the object nz_height) onto the map_nz2 object with tm_dots() function.
4.1.1 Example 1: What state contains the geographic centroid of tornado activity?
Import the tornado data. Remove Alaska, Hawaii, and Puerto Rico tornadoes. The native CRS is latitude/longitude so we transform it to a Web Mercator (EPSG 3857) used by Google Maps, OpenStreetMap, Bing, ArcGIS, ESRI.
Alltors.sfdf <- read_sf(dsn = "1950-2018-torn-initpoint") %>%
filter(!st %in% c("HI", "AK", "PR")) %>%
filter(yr >= 1994) %>%
st_transform(crs = 3857)What is the geometry type in the simple feature column (sfc)?
## Geometry set for 30497 features
## geometry type: POINT
## dimension: XY
## bbox: xmin: -13855710 ymin: 2820573 xmax: -7548241 ymax: 6331042
## CRS: EPSG:3857
## First 5 geometries:## POINT (-9625796 3665259)## POINT (-9075878 3214975)## POINT (-9085897 3246452)## POINT (-9176066 3480433)## POINT (-9170500 3493271)The start location of the tornado track has a POINT geometry. The CRS specified by the EPSG number specifies a particular proj4string.
We use the st_combine() function to create a single feature column with a MULTIPOINT geometry. All the genesis locations are geometrically combined into a single MULTIPOINT simple feature (no attributes).
## Geometry set for 1 feature
## geometry type: MULTIPOINT
## dimension: XY
## bbox: xmin: -13855710 ymin: 2820573 xmax: -7548241 ymax: 6331042
## CRS: EPSG:3857## MULTIPOINT ((-9625796 3665259), (-9075878 32149...Next we use the st_centroid() function to find the geographic center of the simple feature column.
The result is a single feature of type POINT. In this case the order of operations is cummutative, but it is not in general.
To create a map, we first get a simple feature data frame (geometry type: MULTIPOLYGON) of the continental U.S. state borders from the spData package.
library(USAboundaries)
States <- us_states() %>%
filter(!state_name %in% c("Alaska", "Hawaii", "Puerto Rico", "District of Columbia"))
st_crs(States)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCS["WGS 84",
## DATUM["WGS_1984",
## SPHEROID["WGS 84",6378137,298.257223563,
## AUTHORITY["EPSG","7030"]],
## AUTHORITY["EPSG","6326"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4326"]]## Geometry set for 48 features
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -124.7258 ymin: 24.49813 xmax: -66.9499 ymax: 49.38436
## CRS: EPSG:4326
## First 5 geometries:## MULTIPOLYGON (((-68.92401 43.88541, -68.87478 4...## MULTIPOLYGON (((-114.7997 32.59362, -114.8094 3...## MULTIPOLYGON (((-94.61792 36.49941, -94.3612 36...## MULTIPOLYGON (((-75.77379 39.7222, -75.75323 39...## MULTIPOLYGON (((-85.60516 34.98468, -85.47434 3...CAUTION: The spData package contains a spatial data frame for getting state-level boundaries with the same name (us_states)!!
We note that the CRS of the state borders does not match the CRS of the tornadoes so we use the st_transform() function with the argument crs = st_crs(Alltors.sfdf).
Which state contains the centroid? The function st_contains() returns a TRUE for the state that contains the center.
The result is a 48 x 1 matrix with all entries FALSE except the state containing the centroid.
Select the geometry containing the centroid.
Then to make a map using the functions from the tmap package we start with the state borders group, then add the tornadoes as a group with a dot layer, then add the state borders subsetted by the state containing the centroid, finally add the centroid location as a group with bubbles as a layer.
tm_shape(States) +
tm_borders() +
tm_shape(Alltors.sfdf) +
tm_dots() +
tm_shape(States[mtrx, 1]) +
tm_borders(col = "red") +
tm_shape(centerState) +
tm_bubbles(size = .3, col = "red")
Note that we need to think about the order of the layers.
4.1.2 Example 2: State-level tornado rates
Here we compute and display the annual tornado occurrence rate (per 10,000 square km) for each state.
Start by first determining the intersections of the state polygons and the tornado points.
## [1] 48 30497The result is a 48 (states) by 30,497 (tornadoes) matrix of logical values indicating whether or not each tornado occurred within each state.
Next we use the rowSums() function to get the total number of TRUE entries for each states.
## [1] 54 101 984 14 831 1054 225 1374 100 1387 1252 636 978 224 346
## [16] 1186 193 208 71 750 3497 462 577 697 54 619 38 1200 1254 217
## [31] 759 5 470 1603 622 1110 2181 45 46 68 60 263 1271 22 45
## [46] 780 400 12The key is that the order of the elements from the rowSums() function matches the order of the states in States.
Also see the aggregate() function (KBDI.Rmd).
So we include the counts as a column in the States simple feature data frame. We include the area of each state. This allows us to compute the rate per 10,000 sq. km. The spatial unit is meter so areas are in sq. meter. To change square meter to square kilometer we multiply by 10 billion (10^10).
nT <- rowSums(mtrx)
StateArea <- st_area(States)
States <- States %>%
mutate(nT,
rate = nT/StateArea/(2018 - 1994 + 1) * 10^10)
head(States)## Simple feature collection with 6 features and 14 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -12781060 ymin: 3550008 xmax: -7452828 ymax: 6340332
## CRS: EPSG:3857
## statefp statens affgeoid geoid stusps name lsad aland
## 1 23 01779787 0400000US23 23 ME Maine 00 79885221885
## 2 04 01779777 0400000US04 04 AZ Arizona 00 294198560125
## 3 05 00068085 0400000US05 05 AR Arkansas 00 134771517596
## 4 10 01779781 0400000US10 10 DE Delaware 00 5047194742
## 5 13 01705317 0400000US13 13 GA Georgia 00 149169848456
## 6 27 00662849 0400000US27 27 MN Minnesota 00 206232257655
## awater state_name state_abbr jurisdiction_type
## 1 11748755195 Maine ME state
## 2 1027346486 Arizona AZ state
## 3 2960191698 Arkansas AR state
## 4 1398720828 Delaware DE state
## 5 4741100880 Georgia GA state
## 6 18929176411 Minnesota MN state
## geometry nT rate
## 1 MULTIPOLYGON (((-7672586 54... 54 0.12242390 [1/m^2]
## 2 MULTIPOLYGON (((-12779442 3... 101 0.09306545 [1/m^2]
## 3 MULTIPOLYGON (((-10532819 4... 984 1.91609715 [1/m^2]
## 4 MULTIPOLYGON (((-8435099 48... 14 0.63421637 [1/m^2]
## 5 MULTIPOLYGON (((-9529523 41... 831 1.53764865 [1/m^2]
## 6 MULTIPOLYGON (((-10823487 6... 1054 0.91800271 [1/m^2]Then we make a map. Here since all the information is in a single spatial data frame States we only have one group. The group has two layers; borders and fill. The fill layer has the color aesthetic pointed to the column labeled ‘rate’.
tm_shape(States) +
tm_borders(col = "gray70") +
tm_fill(col = "rate",
title = "Annual Rate\n[/10,000 sq. km]") +
tm_layout(legend.outside = TRUE)
Note: functions from tmap (1) expect data as spatial objects rather than data frames and (2) variable names need to be surrounded by quotes.
All functions start with tm_. The first function in a group needs to be tm_shape(). It specifies the spatial object. Functions following tm_shape() specify sequential aesthetics as layers. Layers are divided into base and derived.
Base layers
tm_polygons(): Draws polygons; coltm_symbols(): Draws symbols; size, col, shapetm_lines(): Draws polylines; col, lwdtm_raster(): Draws a raster; coltm_text(): Add text labels; text, size, col
Derived layers
tm_fill(): Fills the polygons; seetm_polygons()tm_borders(): Draws polygon borders; nonetm_bubbles(): Draws bubbles; seetm_symbols()tm_squares(): Draws squares; seetm_symbols()tm_dots(): Draws dots; seetm_symbols()tm_markers(): Draws markers; seetm_symbols()andtm_text()tm_iso()Draws iso/contour lines; seetm_lines()andtm_text()
Each aesthetic (layer) can take a constant value or a data variable name. For instance, tm_fill(col = 'green') colors all polygons green, while tm_fill(col = "var1"), where "var1" is the name of a data variable in the shape object, creates a choropleth. If a vector of constant values or variable names are provided, the output is a set of maps.
The following layers are cartographic elements:
tm_grid(): Add coordinate grid linestm_credits(): Add credits text labeltm_compass(): Add map compasstm_scale_bar(): Add scale bar
For example here we create a choropleth map of countries using an index of happiness.
The simple feature spatial object World is the only required argument in the tm_shape() function. A polygon layer is added where the fill color (col =) is set to the column HPI in the attribute table of the simple feature data frame.
The legend title is set with the title = argument. For example, here we create a title with the expression() and paste() functions then use the title in the tm_fill() function. Here we also use tm_fill() plus tm_borders() rather than the tm_polygons() as the layer.
legend_title <- expression(paste("Area (km", {}^2, ")"))
tm_shape(nz) +
tm_fill(col = "Land_area",
title = legend_title) +
tm_borders()
The default are sensible breaks. The argument breaks = sets the breaks manually. The argument n = sets the number of bins into which numeric variables are categorized. The palette = argument defines the color scheme, for example BuGn from the RColorBrewer package.
4.1.3 Map layout, facets, and inserts
Layout functions help create a cartographic map. Elements include the title, the scale bar, margins, aspect ratios, etc. For example, here elements such as a north arrow and a scale bar are added with tm_compass() and tm_scale_bar(), respectively and the tm_layout() function is used to add the title and background color.
map_nz +
tm_compass(type = "8star",
position = c("left", "top")) +
tm_scale_bar(breaks = c(0, 100, 200),
text.size = 1) +
tm_layout(title = "New Zealand",
bg.color = "lightblue")
Faceted maps (referred to as ‘small multiples’) are composed of several maps arranged side-by-side. Facets enable the visualization of how spatial relationships change with respect to another variable.
For example, here the faceted variable is time (year). The simple feature data frame is from the spData package. We first filter the data frame keeping only the years 1970, 1990, 2010, and 2030.
Note: The operator %in% acts like a recursive or. If year == 1970 or year == 1990, … For example,
## [1] 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983
## [16] 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998
## [31] 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013
## [46] 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028
## [61] 2029 2030 2031## [1] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [37] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [61] FALSE TRUE FALSEWe then make a map.
tm_shape(world) +
tm_polygons() +
tm_shape(urb_1970_2030) +
tm_symbols(col = "black",
border.col = "white",
size = "population_millions") +
tm_facets(by = "year",
nrow = 4,
free.coords = FALSE)The above code chunk demonstrates key features of faceted maps created with functions from the tmap package.
- Shapes that do not have a facet variable are repeated (the countries in world in this case).
- The
by =argument which varies depending on a variable (yearin this case). - nrow/ncol setting specifying the number of rows (and columns) that facets should be arranged into.
- The
free.coords =argument specifies whether each map has its own bounding box.
Small multiples are also generated by assigning more than one value to one of the aesthetic arguments. For example here we map the happiness index (HPI) and gross domestic product per capita (gdp_cap_est).
tm_shape(World) +
tm_polygons(c("HPI", "gdp_cap_est"),
style = c("pretty", "kmeans"),
palette = list("RdYlGn", "Purples"),
title = c("Happy Planet Index", "GDP per capita")) 
Two maps are created each with a different scale. All arguments of the layer functions can be vectorized, one for each small multiple map. Arguments that normally take a vector, such as palette =, are placed in a list().
Multiple map objects can also be arranged in a single plot with the tmap_arrange() function.
4.1.4 Example 3: Tornado locations by month
Suppose we want a map showing the location of tornadoes by month. First add the world map with countries as a polygon layer (fill is gray by default). Then overlay the U.S. country polygons and color them white. Then add the tornado genesis locations as dots faceted by month. Then add the state polygons as borders. Use a North American Lambert Azimuthal Equal Area Projection and make this the master layer. Finally add cartographic elements.
tm_shape(world) +
tm_polygons() +
tm_shape(world[world$name_long == "United States", ]) +
tm_polygons(col = "white") +
tm_shape(Alltors.sfdf) +
tm_dots() +
tm_facets(by = "mo", as.layers = TRUE) +
tm_shape(States, is.master = TRUE) +
tm_borders() +
tm_compass(type = "arrow", position = c("left", "bottom")) +
tm_scale_bar(position = c("left", "bottom"), text.size = .75) +
tm_style("natural") +
tm_layout(main.title = "Contiguous U.S. Tornadoes [1950-2018]",
main.title.position = "center", main.title.size = .85,
panel.labels = month.name) +
tm_credits(c(rep("", 11), "Data Source: U.S. SPC"),
position = c("right", "bottom"))We can clearly see the spread of tornado activity from the southeast in January to the central Plains in May to the northern Plains during July and August and then back to the southeast in December.
4.1.5 Inset map
An inset map is a smaller map set within or next to the main map. An inset map is used to contextualize the geographic study area. Here we create a map of the central part of New Zealand’s Southern Alps. The inset map shows where the main map is in relation to all of New Zealand.
The first step is to define the area of interest. Here it is done here by creating a new spatial object nz_region using the st_bbox() function and the st_as_sfc() to make it a simple feature column.
nz_region <- st_bbox(c(xmin = 1340000, xmax = 1450000,
ymin = 5130000, ymax = 5210000),
crs = st_crs(nz_height)) %>%
st_as_sfc()The second step is to create a base map showing New Zealand’s Southern Alps area. This is the closeup view of where the most important message is stated. The region is clipped to the simple feature column nz_region created above. The layers include a raster of elevations and locations of high points. A scale bar is included.
nz_height_map <- tm_shape(nz_elev,
bbox = nz_region) +
tm_raster(style = "cont",
palette = "YlGn",
legend.show = TRUE) +
tm_shape(nz_height) +
tm_symbols(shape = 2,
col = "red",
size = 1) +
tm_scale_bar(position = c("left", "bottom"))
nz_height_map
The third step is to create the inset map. It gives a context and helps to locate the area of interest. This map clearly indicates the location of the main map.
nz_map <- tm_shape(nz) +
tm_polygons() +
tm_shape(nz_height) +
tm_symbols(shape = 2,
col = "red",
size = .1) +
tm_shape(nz_region) +
tm_borders(lwd = 3)
nz_map
The final step is to combine the two maps. The viewport() function from the grid package is used to give a center location (x and y) and the size (width and height) of the inset map.
Additional details are available here: https://geocompr.robinlovelace.net/adv-map.html
4.1.6 Turn a static map into an interactive map
Leaflet is a popular open-source platform (javascript) for making/serving interactive maps. The leaflet package provides access to the platform.
As an example, here we create a world map by using the default mapping arguments in leaflet() and zoomable tiles with the addTiles() layer. The piping operator (%>%) is used to ‘add’ layers.
Note: if the map doesn’t render in your Viewer tab, then click the Viewer tab and select “Show in new window”.
Here we use setView() with longitude and latitude arguments in decimal degrees to center the map on FSU and use a zoom of 17.
A cool feature of tmap is that we can create an interactive map using the same code as we used to create a static map.
For example our static map of New Zealand (map_nz) is viewed interactively by switching to view mode.
## tmap mode set to interactive viewingWith the interactive mode turned on, all maps produced with tmap will launch as zoomable HTML. This feature includes the ability to specify the basemap with tm_basemap() (or tmap_options()) as demonstrated here.
Q: Why is there no topography for New Zealand?
We can also create interactive maps with the tmap_leaflet() function.
The view mode in tmap works with faceted plots. The argument sync in tm_facets() is used to produce multiple maps with synchronized zoom and pan settings.
world_coffee <- left_join(world,
coffee_data,
by = "name_long")
tm_shape(world_coffee) +
tm_polygons(c("coffee_production_2016",
"coffee_production_2017")) +
tm_facets(nrow = 1, sync = TRUE)Change the view mode back to plot.
## tmap mode set to plotting4.1.7 Example 4: Tornado occurrences on a grid
Import the tornado data convert it to geographic coordinates.
Alltors2.sfdf <- st_read(dsn = "1950-2018-torn-aspath") %>%
filter(!st %in% c("HI", "AK", "PR")) %>%
filter(yr >= 2014 & mag >= 0)## Reading layer `1950-2018-torn-aspath' from data source `/Users/jelsner/Desktop/Projects/ASSFG-Book/1950-2018-torn-aspath' using driver `ESRI Shapefile'
## Simple feature collection with 63645 features and 22 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: -163.53 ymin: 18.13 xmax: -64.9 ymax: 61.02
## CRS: 4326Create the leaflet map with a view set on campus. First use the addTiles() function and then the addPolylines() function with the argument data = Alltors2.sfdf. The piping operator is used to add layers.
m <- leaflet() %>%
setView(-84.29849, 30.44188, zoom = 17) %>%
addTiles() %>%
addPolylines(data = Alltors2.sfdf)
mWe can also use the mapView() function from the mapview package. It also provides interactivity for easy and quick visualization during spatial data analysis. But it is not intended for fine-tuned presentation quality map production.